package com.shsmi.util;
//球面(经纬度、大地)坐标 与 平面坐标  互转
//1:100万图幅

import com.shsmi.entity.CityPoint;
import com.shsmi.entity.GeoPoint;

public class GaussXYDeal {

    public static GeoPoint GaussToBL(double X, double Y)   //, double *longitude, double *latitude)
    {
        int ProjNo;
        int ZoneWide;           //带宽
        double[] output = new double[2];
        double latitude0, longitude1, latitude1, longitude0, X0, Y0, xval, yval;//latitude0,
        double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
        iPI = 0.0174532925199433;          // 3.1415926535898/180.0;
        a = 6378137.0;



        f = 1.0 / 298.257223563;    //54年北京坐标系参数
//        a = 6378140.0;
//        f = 1 / 298.257;           //80年西安坐标系参数
//        a = 6378137.0;
//        f = 1 / 298.257222101; //2000椭球参数


        ZoneWide = 3;                       //6度带宽
        ProjNo = (int) (X / 1000000L);        //查找带号
//        longitude0 = (ProjNo - 1) * ZoneWide + ZoneWide / 2;

        //longitude0 = 121.4669444;
        longitude0 = 121.467168475;
        longitude0 = longitude0 * iPI;     //中央经线

//        latitude0 = 31.23527778;
        latitude0 = 31.2354333222222;
        latitude0 = latitude0 * iPI;     //中央经线

//        X0 = ProjNo * 1000000L + 500000L;
        X0 = 0;
        Y0 = 0;
        xval = X - X0;
        yval = Y - Y0;           //带内大地坐标
        e2 = 2 * f - f * f;
        e1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2));
        ee = e2 / (1 - e2);
        M = yval;
        u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
        fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.sin(
                4 * u)
                + (151 * e1 * e1 * e1 / 96) * Math.sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.sin(8 * u);
        C = ee * Math.cos(fai) * Math.cos(fai);
        T = Math.tan(fai) * Math.tan(fai);
        NN = a / Math.sqrt(1.0 - e2 * Math.sin(fai) * Math.sin(fai));
        R = a * (1 - e2) / Math.sqrt((1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin
                (fai) * Math.sin(fai)));
        D = xval / NN;
        //计算经度(Longitude) 纬度(Latitude)  弧度
        longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
                * D * D * D * D / 120) / Math.cos(fai);
        latitude1 = latitude0 + fai - (NN * Math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
                + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
        //转换为度 DD
        output[0] = longitude1 / iPI;
        output[1] = latitude1 / iPI;
        GeoPoint geoPoint = new GeoPoint();
        geoPoint.setLat(output[0]);
        geoPoint.setLon(output[1]);
        return geoPoint;
    }

    //由经纬度反算成高斯投影坐标
    public static double[] BLToGauss2(double longitude, double latitude) {
        int ProjNo = 0;
        int ZoneWide;             //带宽
        double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
        double a, f, e2, ee, NN, T, C, A, M, iPI;
        iPI = 0.0174532925199433;               //3.1415926535898/180.0;
        ZoneWide = 3;                           //6度带宽
//        a = 6378245.0;
//        f = 1.0 / 298.3;               //54年北京坐标系参数
//        a = 6378140.0;
//        f = 1 / 298.257; //80年西安坐标系参数
        a = 6378137.0;
        f = 1 / 298.257222101; //2000椭球参数
        ProjNo = (int) (longitude / ZoneWide);
//        longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
        longitude0 = 121.4669444;
        longitude0 = longitude0 * iPI;
        latitude0 = 31.23527778;
        latitude0 = latitude0 * iPI;
        System.out.println(latitude0);
        longitude1 = longitude * iPI; //经度转换为弧度
        latitude1 = latitude * iPI; //纬度转换为弧度
        e2 = 2 * f - f * f;
        ee = e2 * (1.0 - e2);
        NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1 - latitude0) * Math.sin(latitude1 - latitude0));
        T = Math.tan(latitude1 - latitude0) * Math.tan(latitude1 - latitude0);
        C = ee * Math.cos(latitude1 - latitude0) * Math.cos(latitude1 - latitude0);
        A = (longitude1 - longitude0) * Math.cos(latitude1 - latitude0);
        M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * (latitude1 - latitude0) - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2
                * e2 / 1024) * Math.sin(2 * (latitude1 - latitude0))
                + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * (latitude1 - latitude0)) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * (latitude1 - latitude0)));
        xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
        yval = M + NN * Math.tan(latitude1 - latitude0) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
                + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
        X0 = 1000000L * (ProjNo + 1) + 500000L;
//        X0 = 0;
        Y0 = 0;

        double[] lonlat = new double[2];
        lonlat[0] = xval + X0;
        lonlat[1] = yval + Y0;
        return lonlat;
    }

    //变更后的由经纬度反算成高斯投影坐标, meridian 为mesh中央经度
    public static CityPoint BLToGauss(double longitude, double latitude) {
//        double meridian = 121.4644444;//中央经线
//        double meridian = 121.4669444;//中央经线
//
//        double zeroLon = 121.4669444;//原点经度
//        double zeroLat = 31.23527778;//原点纬度

        double meridian = 121.467168475;//中央经线

        double zeroLon = 121.467168475;//原点经度
        double zeroLat = 31.2354333222222;//原点纬度

        int ProjNo = 0;
        double longitude1, latitude1, X0, Y0, xval, yval;
        double a, f, e2, ee, NN, T, C, A, M, iPI;
        iPI = 0.0174532925199433; //3.1415926535898/180.0;
        a = 6378137.0;
        f = 1.0 / 298.257223563; //54年北京坐标系参数
//        a = 6378140.0;
//        f = 1 / 298.257; //80年西安坐标系参数
//        a = 6378137.0;
//        f = 1 / 298.257222101; //2000椭球参数
        meridian = meridian * iPI;
        longitude1 = longitude * iPI; //经度转换为弧度
        latitude1 = latitude * iPI; //纬度转换为弧度
        zeroLon = zeroLon * iPI;
        zeroLat = zeroLat * iPI;

        e2 = 2 * f - f * f;
        ee = e2 * (1.0 - e2);
        NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1 - zeroLat) * Math.sin(latitude1 - zeroLat));
        T = Math.tan(latitude1 - zeroLat) * Math.tan(latitude1 - zeroLat);
        C = ee * Math.cos(latitude1 - zeroLat) * Math.cos(latitude1 - zeroLat);
        A = (longitude1 - meridian) * Math.cos(latitude1 - zeroLat);
        M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * (latitude1 - zeroLat) - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2
                * e2 / 1024) * Math.sin(2 * (latitude1 - zeroLat))
                + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * (latitude1 - zeroLat)) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * (latitude1 - zeroLat)));
        xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
        yval = M + NN * Math.tan(latitude1 - zeroLat) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
                + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
//        X0 = 1000000L * (ProjNo + 1) + 500000L;
        X0 = 0;//修改
        Y0 = 0;
        double[] lonlat = new double[2];
        lonlat[0] = xval + X0;
        lonlat[1] = yval + Y0;
        CityPoint cityPoint = new CityPoint();
        cityPoint.setX(lonlat[0]);
        cityPoint.setY(lonlat[1]);
        return cityPoint;
    }

    public static void main(String[] args) throws Exception {
        GaussXYDeal g = new GaussXYDeal();
//        double[] xy = g.BLToGauss(121.4669444, 31.23527778);
//        System.out.println(xy[0]);
//        System.out.println(xy[1]);
//        System.out.println();

//        CityPoint cityPoint = g.BLToGauss(121.467155, 31.23545044);
        CityPoint cityPoint = g.BLToGauss(121.719839, 30.97348655);

        System.out.println("x: " + cityPoint.getX());
        System.out.println("y: " + cityPoint.getY());

        GeoPoint geoPoint = g.GaussToBL(cityPoint.getX(), cityPoint.getY());
        System.out.println(geoPoint.getLat());
        System.out.println(geoPoint.getLon());

        System.out.println();

 
        /*double[] lonlat = g.GaussToBL(xy[0], xy[1]);
        System.out.println("lon: "+lonlat[0]);
        System.out.println("lat: "+lonlat[1]);
 
        System.out.println();
        double[] a = g.GaussToBLToGauss(113.26846,23.55979, 113.3125);
        System.out.println("x: "+a[0]);
        System.out.println("y: "+a[1]);*/
    }

}